

% Define useful functions

bargaining = Params.barg ;

Sbar  = @(zeta) Params.tau ./ ( zeta - 1 ) ./ ( Params.tau + zeta ) ;

SbarPrimeOverSbar = ...
        @(zeta) - ( Params.tau + 2 * zeta - 1 ) ./ ( zeta - 1) ./ ( Params.tau + zeta );

Phi   = @(v,zeta) Params.rho / bargaining * v ./ ( ( Params.b + v ) .* Sbar(zeta) ) ;
PhiVOverPhi  = @(v) Params.b ./ v ./ ( Params.b + v ) ;
PhiZOverPhi  = @(zeta) - SbarPrimeOverSbar(zeta) ;

Theta = @(v,zeta) ( 1 ./ Params.m0 .* Phi(v,zeta) .* ( ( Params.b + v ) / Params.B0 ).^zeta ).^( 1 / ( 1 - Params.alpha ) )  ;

uRate0 = @(v,zeta) ( Params.D + Params.delta * zeta ) ./ ( Params.D + Params.delta .* zeta + Phi(v,zeta) ) ;
uRate  = @(v,zeta) min( max( uRate0(v,zeta) , Grids.minu ) , Grids.maxu ) ;
uV = @(v,zeta) uRate(v,zeta) ./ ( Params.DLearn + Params.delta * zeta + Phi(v,zeta) ) ...
               .* Phi(v,zeta) .* PhiVOverPhi(v) ;
uZZ = @(v,zeta) uRate(v,zeta) ./ ( Params.DLearn + Params.delta * zeta ) ...
               .* ( Params.delta * ( 1 - uRate(v,zeta) ) + Phi(v,zeta) .* PhiZOverPhi(zeta) ) ;

uMean = @(v,zeta,L) sum( L .* Grids.mx .* uRate(v,zeta) ) / sum( L .* Grids.mx ) ;
           
InstantSep = ...
        @(v,zeta) ( Params.B0 ./ ( Params.b + v ) ).^zeta ;

Phi2  = @(v,zeta) Params.rho / bargaining * v ./ ( ( Params.b + v ) .* Sbar(zeta) ) ...
               .* InstantSep(v,zeta) ;
                       
Stilde = ...
        @(zeta) Params.tau * ( Params.tau * ( zeta + Params.kappa - 1 ) + zeta.^2 + ( Params.kappa - 1 ) * ( zeta + Params.kappa) ) ...
             ./ (  ( Params.kappa -1) * ( Params.tau + Params.kappa ) * ( zeta - 1) .* ( Params.tau + zeta )  ) ;
                        
ExpProd = ...
        @(v,zeta) Params.zl0 ./ Params.rho ...
               .* Params.kappa .* zeta ./ ( Params.kappa - 1 ) ./ ( zeta - 1 ) ;
EpZOverEp = @(zeta) 1 ./ zeta ./ ( zeta - 1 ) ;         
ExpProdZ  = @(v,zeta) ExpProd(v,zeta) .* EpZOverEp(zeta) ;
           
         
G = ...
    @(v,zeta) Params.omega * uRate(v,zeta) .* Params.b ...
            + Params.omega * ( 1 - uRate(v,zeta) ) .* ( Params.b + v ) ...
              .* ( 1 - bargaining + bargaining .* ExpProd(v,zeta) ) ...
            + Params.psi * ( 1 - uRate(v,zeta) ) .* ( Params.b + v ) ...
              .* ExpProd(v,zeta) ;
          
LocalHC = @(v,zeta) ...
    ( Params.DLearn ./ ( Params.DLearn + Params.phi * uRate(v,zeta) ) ).^( ( Params.omega + Params.eps * ( 1 + Params.eta ) ) / Params.P1 ) ;

% define two different profit functions: one for the return to entry,
% and one for vacancy creation. see Online Appendix for computation
Profits = ...
    @(v,zeta,x)    (    Params.B0.^zeta .* ( Params.b + v ).^( 1 - zeta ) ...
                     .* v.^( - Params.alpha ) .* Sbar(zeta) ...
                   ).^( 1 / ( 1 - Params.alpha ) ) ...
                .* x .* LocalHC(v,zeta) ...
                .* (    ( Params.b + v ) ...
                     .* G(v,zeta).^Params.eps ...
                   ).^( - Params.psi / Params.P1 ) ;
% \hat{J} in paper
               
ProfitsVac = @(v,zeta,x,tax) Profits(v,zeta,x) .* tax.^Params.tax ;
% \tilde{J} in paper

Vacancies = ...
    @(v,zeta,x,tax) ProfitsVac(v,zeta,x,tax).^Params.gamma ;
% use \tilde{J} as in paper

AdjxFunc = ...
         @(x) exp( Params.stdaCondx^2 / 2 / ( Params.omega + Params.eps * ( 1 + Params.eta ) )^2 ) ...
         * x.^( Params.aCondx / ( Params.omega + Params.eps * ( 1 + Params.eta ) ) ) ;
              
              
              
GV = @(v,zeta)   Params.omega * Params.b * uV(v,zeta)  ...
               + Params.omega * ( 1 - uRate(v,zeta) )          .* ( 1 - Params.barg + Params.barg * ExpProd(v,zeta) ) ...
               - Params.omega * ( Params.b + v ) .* uV(v,zeta) .* ( 1 - Params.barg + Params.barg * ExpProd(v,zeta) ) ...
               + Params.psi   * ( 1 - uRate(v,zeta) )          .* ExpProd(v,zeta) ...
               - Params.psi   * ( Params.b + v ) .* uV(v,zeta) .* ExpProd(v,zeta) ;

GZ = @(v,zeta)   Params.omega * Params.b * uZZ(v,zeta)  ...
               + Params.omega * ( 1 - uRate(v,zeta) )  .* ( Params.b + v ) .* Params.barg .* ExpProdZ(v,zeta) ...
               - Params.omega * uZZ(v,zeta)            .* ( Params.b + v ) .* ( 1 - Params.barg + Params.barg * ExpProd(v,zeta) ) ...
               + Params.psi   * ( 1 - uRate(v,zeta) ) .* ( Params.b + v ) .* Params.barg .* ExpProdZ(v,zeta) ...
               - Params.psi   * ( Params.b + v ) .* uZZ(v,zeta) .* ExpProd(v,zeta) ;
           
           
Lhat              = ...
    @(v,zeta,x)    ( Params.b + v ).^( ( 1 + Params.eta + Params.psi ) / Params.P1 ) ...
                .* x.^( ( 1 + Params.eta - Params.omega ) / ( Params.omega + Params.eps * ( 1 + Params.eta ) ) ) ...
                ./ G(v,zeta).^( ( Params.omega + Params.psi ) / Params.P1 ) ...
                .* AdjxFunc(x) ;
                
UnitProfits       =  @(v,zeta,x) ...
                   Profits(v,zeta,x) ...
                .* Lhat( v , zeta , x ) .* uRate( v , zeta ) .* Theta( v , zeta ) ;
% use \hat{J} (only used for returns to entry)

MeanHC = @(x0,Lh,v,zeta) ...
    (   ( 1 - x0 ) ...
      / ( 1 - x0 * sum( Lh ./ ( 1 + Params.phi / Params.DLearn * uRate(v,zeta) ) .* Grids.mx ) ...
                 / sum( Lh .* Grids.mx ) ) ...
     ).^Params.P2 ;
    

            
K0               = @(zeta) 1 - bargaining + ...
                            bargaining * Params.zl0 ...
                          * Params.kappa / ( Params.kappa - 1 ) ...
                          * zeta ./ ( zeta - 1 ) ;
                       
                        
                        
 Es = @(zeta) ...
     Params.tau * ( zeta.^2 + zeta .* ( Params.kappa + Params.tau - 1 ) ...
                    + ( Params.kappa - 1 ) * ( Params.kappa + Params.tau ) ...
                   ) ...
   ./ ( ( zeta - 1 ) .* ( Params.kappa - 1 ) ...
        .* ( zeta + Params.tau ) .* ( Params.kappa + Params.tau ) ) ;

 
NomWagesNoHC = @(v,zeta) ...
       Grids.x .* ( Params.b + v ).^( 1 - Params.psi / Params.P1 ) ...
    .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) .* ( 1 - bargaining + bargaining * ExpProd(v,zeta) ) ;            

NomWagesHC = @(v,zeta) ...
       Grids.x .* ( Params.b + v ).^( 1 - Params.psi / Params.P1 ) ...
    .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) .* ( 1 - bargaining + bargaining * ExpProd(v,zeta) ) ...
    .* LocalHC(v,zeta) ;            

              
% to update entry condition
EntryRatio = @(x0,v,zeta,Grids,barg,GEC) ...
           sum( UnitProfits(v,zeta,Grids.x) .* Grids.mx ) ...
           * sum( Lhat(v,zeta,Grids.x) .* Grids.mx ).^( (1+Params.gamma)*Params.psi/(1+Params.eta+Params.psi) ) ...
           * MeanHC(x0,Lhat(v,zeta,Grids.x),v,zeta)^(1+Params.gamma) ...
           / ( GEC * Params.ce * ( (1-barg)*barg^(Params.alpha/(1-Params.alpha)) )^(-1-Params.gamma) ) ;

% local production, income, profits and land rents
HomeProdHC = @(v,zeta) ...
       Grids.x .* ( Params.b + v ).^( 1 - Params.psi / Params.P1 ) ...
    .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) .* Params.b ...
    .* LocalHC(v,zeta) ;    

FlowProfitsHC = @(v,zeta) ...
       Grids.x .* ( Params.b + v ).^( 1 - Params.psi / Params.P1 ) ...
    .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) .* ( 1 - bargaining ) .* ( ExpProd(v,zeta) - 1 ) ...
    .* LocalHC(v,zeta) ;

LocalOutputHC = @(v,zeta,Grids) ...
    Grids.x .* ( Params.b + v ).^( 1 - Params.psi / Params.P1 ) ...
    .* G(v,zeta).^( 1 - Params.eps * Params.psi / Params.P1 ) .* ExpProd(v,zeta) ...
    .* LocalHC(v,zeta) ;    

% aggregate production, income and profits
AgLabIncome = @(v,zeta) ...
    sum( (   uRate(v,zeta)         .* HomeProdHC(v,zeta) ...
           + ( 1 - uRate(v,zeta) ) .* NomWagesHC(v,zeta) ...
         ) .* Lhat(v,zeta,Grids.x) .* Grids.mx ) ...
  / sum( Lhat(v,zeta,Grids.x) .* Grids.mx ) ;

AgProfits = @(v,zeta,Grids) ...
    sum(    FlowProfitsHC(v,zeta) .* ( 1 - uRate(v,zeta) ) .* Lhat(v,zeta,Grids.x) .* Grids.mx ) ...
  / sum( Lhat(v,zeta,Grids.x) .* Grids.mx ) ;

Tau0 = @(v,zeta) ...
    ( sum( UnitProfits(v,zeta,Grids.x) .*                Grids.mx ) ...
         / sum( UnitProfits(v,zeta,Grids.x) .* Grids.tax .* Grids.mx ) ...
    )^(1/(1+Params.gamma)) ;
   
AgTaxExp = @(tau0,phi,v,zeta,Grids) ...
    sum(    FlowProfitsHC(v,zeta) .* ( tau0 * phi - 1 ) ...
         .* ( 1 - uRate(v,zeta) ) .* Lhat(v,zeta,Grids.x) .* Grids.mx ) ...
  / sum( Lhat(v,zeta,Grids.x) .* Grids.mx ) ;

AgOutput = @(v,zeta,Grids) ...
    sum( ( 1 - uRate(v,zeta) ) .* LocalOutputHC(v,zeta,Grids) ...
         .* Lhat(v,zeta,Grids.x) .* Grids.mx ) ...
  / sum( Lhat(v,zeta,Grids.x) .* Grids.mx ) ;

                        
                        